"""
========================================================================
FILE: basis.py

LAST MODIFIED: 20 August 2009

DESCRIPTION: 
Basis classes and constructor for creating and returning a basis 
instance of the required dimensionality and basis
========================================================================
"""
"""
BEGIN LICENSE BLOCK 
Version: MPL 1.1/GPL 2.0/LGPL 2.1

The contents of this file are subject to the Mozilla Public License Version
1.1 (the "License"); you may not use this file except in compliance with
the License. You may obtain a copy of the License at
http://www.mozilla.org/MPL/

Software distributed under the License is distributed on an "AS IS" basis,
WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
for the specific language governing rights and limitations under the
License.

The Original Code is fieldwork.

The Initial Developer of the Original Code is Ju Zhang.
Portions created by the Initial Developer are Copyright (C) 2009
the Initial Developer. All Rights Reserved.

Contributor(s):

Alternatively, the contents of this file may be used under the terms of
either the GNU General Public License Version 2 or later (the "GPL"), or
the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
in which case the provisions of the GPL or the LGPL are applicable instead
of those above. If you wish to allow use of your version of this file only
under the terms of either the GPL or the LGPL, and not to allow others to
use your version of this file under the terms of the MPL, indicate your
decision by deleting the provisions above and replace them with the notice
and other provisions required by the GPL or the LGPL. If you do not delete
the provisions above, a recipient may use your version of this file under
the terms of any one of the MPL, the GPL or the LGPL.

END LICENSE BLOCK
"""

import scipy

# basis matrices of implemented bases for tensor product basis 
basis_matrices_map = {'linear_lagrange':    scipy.array( [[1.0, -1.0],\
												         [0.0,  1.0]] ),\
					 'quadratic_lagrange': scipy.array( [[1.0, -3.0,  2.0],\
														 [0.0,  4.0, -4.0],\
														 [0.0, -1.0,  2.0]] ),\
					 'cubic_hermite':      scipy.array( [[1.0, 0.0, -3.0,  2.0],\
														 [0.0, 1.0, -2.0,  1.0],\
														 [0.0, 0.0,  3.0, -2.0],\
														 [0.0, 0.0, -1.0,  1.0]] ) }

#=======================================================================#

def make_basis( dimensions, basis_types ):
	""" Constructor function for making a basis object 
	basis_types should be a list or tuple of strings
	"""
	
	global implemented_bases

	if basis_types[0] == 'triangular_linear_2d':
		B = triangular_linear_2d()
	elif basis_types[0] == 'triangular_quadratic_2d':
		B = triangular_quadratic_2d()
	else:
		# check for right number of basis_types for specified dimensionality
		if len( basis_types ) != dimensions:
			print 'ERROR: wrong number of basis_types for dimensions'
			return 0
		else:
			basis_matrices = []
			final_type = ''
			for type in basis_types:
				if type not in basis_matrices_map.keys():
					print 'ERROR: unrecognised basis', type
					return 0
				else:
					basis_matrices.append( basis_matrices_map[type] )
					final_type += ' '+type
		
		# instantiate basis object
		B = tensor_product_basis( dimensions, basis_matrices, type )
	
	return B
	
#======================================================================#

class tensor_product_basis:
	""" Tensor product basis object for quad elements
	dimension: number of dimensions
	basis_matrices: a list of tensor product basis matrices, one for each dimension
	type: a string containing the names of basis in each dimension		
	"""
	
	def __init__( self, dimensions, basis_matrices, type ):
		
		self.dimensions = dimensions			# integer of number of dimensions
		self.basis_matrices = basis_matrices	# list of rank 2 arrays
		self.type = type						# string of basis type?
		self.basis_orders = [ (b.shape[0] - 1) for b in self.basis_matrices ]
	
	def eval( self, x ):
		# evaluates the basis function at the given element coordinates
		# returns a list of values to be dot multiplied with nodal values
		# x should be a list of lists of xi coordinates
		# e.g. [ xi1 , xi2 ,...]

		# check x to be of right dimensionality
		if len( x ) != self.dimensions:
			print 'ERROR: wrong number of xi, need', self.dimensions
			return None
		else:		
			# get phis for each dimension
			phid = []
			for i in range( self.dimensions ):
				x_vector = [ x[i]**p for p in range( self.basis_orders[i] + 1 ) ]
				phid.append( scipy.dot( self.basis_matrices[i], x_vector ) )
			
			# do tensor product using kron
			# reverse dimensions to get same implementation as cmiss
			phi = phid[0]
			if self.dimensions > 1:
				for d in range( 1, self.dimensions ):
					phi = scipy.kron( phid[d], phi )
			
			return phi

#======================================================================#

class triangular_linear_2d:
	""" Triangular element bases are hard coded for now
	"""
	
	def __init__( self ):
		self.dimensions = 2
		self.type = 'triangular_linear_2d'
		self.basis_order = 1

	def eval( self, x ):
		phi = [0.0]*3
		phi[0] = 1 - x[0] - x[1]
		phi[1] = x[0]
		phi[2] = x[1]
		
		return phi

class triangular_quadratic_2d:
	""" Triangular element bases are hard coded for now
	"""
	
	def __init__( self ):
		self.dimensions = 2
		self.type = 'triangular_quadratic_2d'
		self.basis_order = 2

	def eval( self, x ):
		phi = [0.0]*6
		phi[0] = ( 1- x[0] - x[1] ) * ( 2 * ( 1 - x[0] - x[1] ) - 1 )
		phi[1] = 4 * ( 1 - x[0] - x[1] ) * x[0]
		phi[2] = x[0] * ( 2*x[0] - 1 )
		phi[3] = 4 * x[0] * x[1]
		phi[4] = x[1] * ( 2*x[1] - 1 )
		phi[5] = 4 * x[1] * ( 1 - x[0] - x[1] )
		
		return phi

